Mon. Not. R. Astron. Soc. 000, ITHHl (20121 Printed 19 September 2012 (MN JAT^X style file v2.2) 



An Alternative Numerical Method for the Stationary 
Pulsar Magnetospheres in the Force-Free System 

Yohsuke Takamori 1 *, Hirotada Okawa 2 f, Makoto Takamoto 3 J, and Yudai Suwa 4 

1 Osaka City University Advanced Mathematical Institute ( OCAMI), 3-3-138 Sugimoto, Sumiyoshi, Osaka 558-8585, Japan 

2 CENTRA, Departamento de Fisica, Instituto Superior Tecnico, Universidade Tecnica de Lisboa - UTL, Av. Rovisco Pais 1, 
1049 Lisboa, Portugal 

3 Max- Planck- Institut far Kernphysik, Saupfercheckweg 1, 69117 Heidelberg, Germany 

4 Yukawa Institute for Theoretical Physics, Kyoto University, Oiwake-cho, Kitashirakawa, Sakyo-ku, Kyoto 606-8502, Japan 



19 September 2012. Pre-print number: OCU-PHYS-374; AP-GR-102 



f 



ABSTRACT 

Stationary pulsar magnetospheres in the force-free system are governed by the pulsar 
equation. Contopoulos, Kazanas, and Fendt (hereafter CKF) numerically solved the 
pulsar equation and obtained a pulsar magnetosphere model called the CKF solution 
that has both closed and open magnetic field lines in f999. The CKF solution is a 
successful solution, but it contains a poloidal current sheet that flows along the last 
open field line and its physics has not been understood well. In this paper, we suggest 
an alternative method to solve the pulsar equation and construct pulsar magnctoshpere 
models without a current sheet introduced by the CKF solution. In our method, the 
pulsar equation is decomposed into Ampere's law and the force-free condition. We 
numerically solve these equations simultaneously with a fixed poloidal current. As a 
result, we obtain a pulsar magnetosphere model without a current sheet, which is 
similar to the CKF solution near the neutron star and has a jet-like structure at a 
distance along the pole. In addition, wc discuss physical properties of the model. 
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1 INTRODUCTION 

Pulsar magnetoshperes play an important role on activities 
of the pulsars. For example, the magnetic field can extract 
the rotational energy and convert to the radiation. Thanks 
to the strong magnetic field, charged particles can be ac- 
celerated and they would become a pulsar wind. In addi- 
tion, because of the acceleration, photons are emitted from 
those charged particles, which becomes a source of X-ray or 
gamma-ray. Thus, in order to understand the pulsar physics, 
it is important to study the magnetic field structures con- 
structed around pulsars. 

Studies of pulsar magnetopheres have started in the 
context of electromagnetostatics in vacuum. The simplest 
model of the pulsar magnetosphere is dipole in which the 
power of r adiation from the pulsa r can be explained as dipole 
radiation (|Pacinil 19671 ). However. lGoldreich fc Julian! (| 19691 ) 
pointed out that plasma should be filled around the neutron 
star because of its rapid rotation with the strong magnetic 
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field. Namely, the pulsar magnetoshpere should be deter- 
mined with motion of plasma considered. 

The system becomes highly non-linear, and therefore 
it is difficult to determine self-consistently both the mo- 
tion of plasma and the configuration of electromagnetic 
field. The simplest assumption is that the kinetic energy 
of plasma is much smaller than the energy of electromag- 
netic field, that is, the inertia of plasma is ignored, which 
is known as the force-free approximation. Because pulsars 
have strong magnetic fields, this assumption would be suf- 
ficiently satisfied. Moreover, when the system is stationary 
and axisymmetric, the electromagnetic field is represented 
by three quantities: the magnetic flux ', the poloidal cur- 
rent 7p(^), and the angular velocity of the magnetic field 
line fiF(vl'). The pulsar magnetosphere in the force-free sys- 
tem is determined by the pu l sar equation (iMichell Il973a ; 
IScharlemann fc Wagoneill 19731 ; lMestellll973l ; I0kamotolll974r ) 
which is a quasi-linear elliptic-type differential equation for 
with /p(^) and ^(v^) as its source terms. For the zero- 
poloidal curr ent, a solution with dipole geometry ha ve been 
constructed (|Mestel fc Wand 1 19791 : IMichell 1 1973b! ) . Some 
cases w ith a non-zero poloidal current have been also dis- 
cussed (iMichejl 1973d l SlBlandfordl 19761 : [Beskin ^ et al.lll983l; 
ISulkanen fc Lovelace! Tl990l : iLvubarskiil ll990l ; iFendt et all 
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ll995l ; lBeskin fc MaIvshkinlll99Sl ). Although they could con- 
struct a pulsar mangetoshpere filled with plasma, due to the 
singularity of the pulsar equation so-called light cylinder, 
some solutions can be applied only inside the light cylinder. 
The others pass through the light cylinder smoothly, but it 
is hard to apply them to realistically astrophysical situation. 

In order to construct realistic pulsar mag- 
netosphre models across the light cylinder, 
IContopoulos. Kazanas. fe Fendtl (1 19991 ) suggested a new 
numerical method which is called the CKF method and 
obtained a solution, the so-called CKF solution which has a 
dipole geometry near the neutron star and passes through 
the light cylinder smoothly. After their work, other proper- 
ties of the CKF solution (e.g., the drift velocity, the e nergy 
loss rate) have been discussed (lOgura fc Koiima 2003 ; 



Goodwin. Mestel. Mestel. fc Wrightl |2004| ; iGruzinovl 12005 ; 
Timokhinl 120061 ). Interestingly, dynamical simulations of 



force-free electrodynamics or magnetohydrodynamics have 
shown that a CKF-like configuration appears as the stead y 
state (|Komissarovll2006l ; [McKinnevi [20061 ; ISpitkovsky||2006h . 
Moreover, the CKF method was extended to the more 
general case that the angular velocity of the magnetic field 
line SIf is not constant and constructed magnetosphere's 
models (|Contopoulosl 120051 ; iTimokhinl 120071 ) . 

The CKF method is a successful and historical mile- 
stone to solve the pulsar equation, and it has been used to 
study pulsar physics. However, there still remains a problem. 
The obtained solution by the CKF method has a poloidal 
current sheet that flows along the equator and the last open 
field line (hereafter, we simply call it "current sheet"). This 
current sheet is ar t ificial ly added to make the current system 
closed. lUzdenskvl (|2003l ) focused on a local region near the 
Y-point, which is the point where the current sheet splits in 
two currents, and imposed that the Maxwell stress at the last 
open field line should be continuous. Consequently, it was 
pointed out that the electromagnetic field becomes infinity if 
the Y-point locates on the light cylinder. In other words, the 
Y-point should locate inside the light cylinder to avoid this 
infinity if one admits the current shee t. Such cases have been 
studied in lGoodwin et all (|2004l ) and ITimokhinl (|2006l ), and 
they obtained solutions that smoothly pass through the light 
cylinder. However, the electromagnetic field has the mini- 
mum energy in th e case that the Y -point locates exactly on 
the light cylinder (|Timokhinll2006l ) so that the solution with 
the minimum energy is unphysical in terms of Uzdensky's 
result. To sum up, it is worth to construct a magnetosphere 
model "without" a current sheet such as introduced in the 
CKF solution. 

Constructing a pulsar magnetosphere model 
without a current sheet h as been studied by 
I Lovelace. Turner, fc Romanoval (|2006l ). Although they 
constructed alternative model of pulsar magnetospheres 
without a current sheet that has a jet flow and a disk wind 
near the equator, they did not construct a CKF-type struc- 
ture w ithout a current sheet. Meanwhile lOgura fc Koiimal 
|2003l) tried to construct a CKF-type structure without a 
current sheet, but they could not obtain such solutions. In 
this paper, we solve the pulsar equation and construct a 
CKF-type structure without a current sheet. It seems that 
the CKF method cannot be used for con structing a CKF- 
type structure without a current sheet (|Ogura fc Koiimal 
120031 ). Therefore, in order to work out our interest, we 



suggest an alternative numerical method to construct pulsar 
magnetosphere models with fixed Ip(9). In addition, we 
study physical properties of pulsar magnetospheres without 
a current sheet. 

This paper is organized as follows. In Sj2j we summarize 
the derivation of the pulsar equation and introduce the CKF 
method. In |3l we explain our method to solve the pulsar 
equation with fixed Ip^) and study pulsar magnetosphere 
models both with and without a current sheet. Then, we 
also study physical properties of our results: the energy loss 
rate; the three-dimensional magnetic field structure in i|4] 
Sj5]is devoted to summary and discussion. 

Through this paperer, we make the speed of light c unity 
and use Gaussian units. 



2 CKF METHOD 

In this section, we briefly summarize the derivation of the 
pulsar equation and the numerical procedure, the so-called 
CKF method, to s olve the pulsar equation suggested by 
IContopoulos et all (| 19991 ). 

Throughout this paper, we focus on stationary and ax- 
isymmetric electromagnetic fields in the force-free system. 
Let E and B be the electric field and the magnetic field, re- 
spectively. Maxwell's equations with stationarity are given 
by 



V • B 

Vx£ 

V ■ E 
V x B 



0. 
0. 

47T/9 e , 

4ttJ, 



(1) 
(2) 
(3) 
(4) 



where p e and J are the electric charge density and the elec- 
tric current density, respectively. In addition, the force-free 
condition is given by 

p c E + JxB = 0. (5) 

Thanks to the symmetries and the force-free condition, the 
electric and the magnetic field can be expressed by the fol- 
lowing three physical quantities; the magnetic flux the an- 
gular velocity of the magnetic field lines f2p, and the poloidal 
electric current Ip. From the force- free condition, we can see 
that f2p and Ip are arbitrary functions of ^ as follows. 

First, let us define the magnetic flux *P and the electric 
current Ip. These are determined by the following surface 
integral; 



ip 



B-dS, 



J-dS, 



(6) 



(7) 



where A is an axisymmetric two dimensional spacelike sur- 
face. Using these functions and the poloidal components of 
Ampere's law, in the the cylindrical coordinates (R,ip,z), 
the magnetic filed B can be written as 



B 



V* x e 



v 2I P (R,z) ^ 



2nR • R (8) 
where e v is a unit vector tangent to the <p direction. Fur- 
thermore, the stationary electric fields can be expressed by 
an electric potential $(-R, 2) as 

£ = V$. (9) 
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Thus, in general, stationary and axisymmetric electromag- 
netic fields can be expressed by ^(R, z), Ip(R,z), and 
&(R, z). In addition, by imposing the force- free condition, 
we find that Ip and fi F become functions of 9. 

Taking the inner product of Eq. ((5| with B, we have 



case, the pulsar equation becomes 



E ■ B = 0. 



(10) 



Substituting Eqs. Q and (J9]) into the above equation, we 
find that the electric potential $ becomes a function of 9, 
that is, $ = Therefore, the electric field can also be 

expressed by \P. 

Next, let us consider the toroidal component of the 
force-free condition. We have 



J P x Bp = 0, 



(11) 



where the subscript P means the poloidal components. From 
the definition of Ip, we have 



Jp 



VJp x e ¥ 
2~kR 



Thus, Eq. |TTJ gives 

{d R Ip)d z <f - (dJ P )d R <f = 0, 



(12) 



(13) 



which implies that Ip is a function of 9. 

Finally, the electric filed E and the magnetic filed B 
can be rewritten in the following form; 



E = 



B 



fi p (tf) 



V*, 



2tt 

V* x e v 2/ P (#) 



2tt7? 



R 



(14) 
(15) 



where we defined that Qf(\I/) := —2'nd<&/d9. Namely, sta- 
tionary and axisymmetric electromagnetic field with the 
force-free condition can be expressed by three scalar func- 
tions *, n F (*), and 1(9). Substituting Eqs. (O and flB) 
into the Gauss's law and Ampere's law, we have 



- / n F (*)V#\ 



-(H) 



= 47tJt, 



(16) 



(17) 



where Jt := J ■ e v . Moreover, we have left the poloidal 
components of the force-free condition. From the poloidal 
components of Eq. ((5}, we obtain the following equation, 



Rtl F p e + 



2I P dip 



Jt- 



(18) 



R 

Combining Eqs. (|16[) , (|17|l . and (|18[) , we find an elliptic type 
differential equation for as follows: 



v. hH+^v*.w- 



i^! Jp ^ =0 (19) 

R 2 1P d y U > l iy i 



where D := 1 — R 2 Q F . This equation is known as the pul- 
sar equation or the Grad-Shafranov equation in the limit 
of strong magnetic field. For pulsar magnetospheres, the 
angular velocity of the magnetic field lines, Q, F , is usually 
assumed as a constant because the magnetic field lines pen- 
etrating the pulsar surface would co-rotate rigidly with the 
pulsar. Here, we also assume that $1 F is a constant. In this 



2RQ F d R 9 + 167r 2 /p^f = 0. 

(20) 

Giving a functional form of Jp(^) and imposing boundary 
conditions, we can obtain a solution of the pulsar equation. 
However, it is well known that a smooth solution of the 
pulsar equation is difficult to obtain due to the existence of 
the light cylinder. 

It is easy to see that the surface where D = 1 — R 2 Q F — 
0, the so-called light cylinder, is singular in the pulsar equa- 
tion. In order to make the derivatives of 9 finite at the light 
cylinder, R = i?LC, we impose one constraint that 



2R LC ti F d R <f\ R=RL 



I 1C 2 r dip 
+ 167T Jp-p- 



= 0. (21) 

If we fix the functional form of Jp($), Eq. (|21|l becomes 
the Neumann boundary condition at the light cylinder, and 
hence the pulsar equation should be solved both inside and 
outside the light cylinder independently. It implies that the 
solutions would not be smooth at the light cylinder in gen- 
eral. 

IContopoulos et al] (|l999h pointed out that Eq. (|21|l can 
be regarded as an equation to determine the functional form 
of Ip(9). In order to obtain a smooth solution beyond the 
light cylinder, they developed the following numerical pro- 
cedure, the so-called CKF method: 

(i) Choose a trial poloidal electric current /p(^). 

(ii) By solving the pulsar equation with the trial Ip both 
inside (— ) and outside (+) the light cylinder, two magnetic 
flux functions, ty^(R, z), are obtained. In general, the value 
of differs from that of at the light cylinder. 

(iii) By substituting $± into Eq. (|21|) . respectively, two 
different poloidal electric currents, Ip±dlp±/ d^ , are ob- 
tained at the light cylinder. Then, a new poloidal electric 
current is determined by 

= mi16^ 2 /p + ^(* + ) + M2l67T 2 / P _^(*-) 

+ #,(*+-*-), (22) 



for 



-(*+ + *-), 



(23) 



where /xi, fJ.2, and ^13 are weight factors satisfying fii+fj,2 = 1 
and ii3 <S 1. 

(iv) Repeat these steps until the difference of becomes 
numerically negligible at the light cylinder. 

Although there is no guarantee whether this scheme con- 
verges, they succeeded in obtaining so-called the CKF so- 
lution. In the CKF method, since Ipdlp/dQ changes in the 
step (iii), the functional form of Ip(9) cannot be determined 
until the calculation finished. Note that one should integrate 
Ipdlp/d^ to obtain Ip(9). At the axis, the electric current 
should satisfy 7p = because there is no electric current 
there. Moreover, Ip also should be zero at the last open 
field line to make the current system closed. However, since 
the equation for Ip is a first order differential equation, we 
can impose one boundary condition. For the CKF solution, 
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the boundary condition as 7p = at the axis is imposed. 
As a result, a current sheet flowing along the last open field 
line is added artificially to make the current system closed. 

One of our purposes is to construct a pulsar magneto- 
sphere model without a current sheet al ong the last open 
field a s introduced in the CKF solution. lOgura fc Koiimal 
(|2003h imposed no current sheet condition in the CKF 
method and carried out their calculation, but they could not 
obtain such solutions. Thus, it seems that the CKF method 
is not useful for obtaining solutions without a current sheet. 
In order that, we want to study the pulsar equation with 
ip(VE') fixed. In the next section, we will suggest an alterna- 
tive numerical method so as to obtain a solution without a 
current sheet and show results. 



3 NUMERICAL STUDY OF PULSAR 

MAGNETOSPHERES WITH FIXED 7 P (*) 

3.1 Numerical procedure 

In order to study the pulsar equation with fixed /p(^I'), let 
us decompose the pulsar equation, (|20[) . into Ampere's law 
and the force-free condition as follows: 



0|* - -d R t> + d 2 * = -S T (R, z), 



(24) 



+ 2RQ 2 F d R ^ - 16vr 2 /pip = -S T (R, z). (25) 

where the prime denotes the derivative with respect to >t and 
we introduced St(R, z) := 8tv 2 R Jt(R, z). The first equation 
is Ampere's law and the second equation is the force-free 
condition with Gauss's law. Both equations can be regarded 
as elliptic type differential equations for the magnetic flux 
Vt(-R, z). Thanks to decomposing the pulsar equation into 
two elliptic type differential equations, we have one more un- 
known function St(R,z) in addition to ^>(R, z) and Ip{^). 
Therefore, we can regard St(R, z) as an adjustable function 
rather than Ip(^). This is the most important point dif- 
ferent from the CKF method, which employed the poloidal 
current Jp(>t) as an adjustable function. In this paper, we 
consider a numerical method for obtaining a solution sat- 
isfying Eqs. Q24p and Q25p simultaneously. Namely, we con- 
sider the problem finding out a set of St) which satisfies 
Eqs. (HU) and (gSJ with fixed 7 P (*). 

There seems to be no singularity like the light cylinder 
in Eqs. (|24|) and (|25[). However, the existence of the light 
cylinder affects the convergence of a numerical iteration. In 
fact, some iterations did not converge if the light cylinder 
existed in the numerical domains. Therefore, to obtain a 
smooth solution across the light cylinder, we construct an 
iterative method as follows: 

(i) Set a function for Ip{^) and give a trial toroidal cur- 
rent <Srtriai(H, z). We fix the functional form of 7p(\E r ) in this 
procedure. 

(ii) By solving Ampere's law (|24[) with the trial 
S'Ttrial {R, z) numerically, we obtain a magnetic flux ^(R, z). 
In general, the obtained ^ and the trial Sr-triai do not satisfy 
the force-free condition (1251) . 

(iii) Using obtained $ and Sr-triai, we make the following 



new trial toroidal current St ucvj (R, z): 
St new (R,z) 

i±t-H^ STin(iM 



+ l-tanh(^) gTout( ^ z) || l e _ D2/(2g2) | 

+ S T L C (I?LC,zK D2/(2t72) , (26) 



where 



S Ti n(R,z) = -{(l + R 2 n 2 F )STtrizl(R,z) 

-2fi|i?aH* + 167r 2 J P Jp}, (27) 

S T out(R,z) = (l + J R 2 «|)- 1 {2S T trial(i?,z) 

+2Q' F Rd R y - 167r 2 / P /p} , (28) 

STLC(RLC,Z) = d 2 R ty\ R = RljC + i?Lc^F9fl*|jJ=ii LC 

- ii?Lcnp 2 (16^ 2 /p/ P )'dfl*| fl=flLC , 

(29) 

and then, 77 and a are constants. We use Smew as the next 

STtrial 

(iv) Repeat these steps until the change of STtraii(i?, z) is 
numerically within a small residual. 

In this routine, we repeat to solve Ampere's law with a 
trial toroidal current that is constructed by using the force- 
free condition. This trial current is given to make this routine 
converge. We performed this routine with some trial currents 
and found that STin (STout) makes this iteration converge if 
the numerical domain is inside (outside) the light cylinder. 
Therefore, by using a hyperbolic tangent and a Gaussian 
function, we make a suitable Smew expected to make this 
iteration converge in a numerical domain including the light 
cylinder. Stiicw almost becomes STin inside the light cylin- 
der, and then the hyperbolic tangent switches Smew from 
STin to STout outside the light cylinder. At the light cylin- 
der, thanks to the Gaussian function, STnew becomes Stlc 
which gives the correct value of St{R, z) at the light cylin- 
der if the light cylinder condition (|2ip is satisfied. We can 
control the behavior of Smew by determining the constants 
a and rj. What values should be adopted to these parameter 
depends on resolution of this calculation. The derivation of 
STnew and the detail are written in AppendixfS] In the step 
(iv), we should determine when we stop this routine. The 
criterion of the convergency is discussed in Appendix [B] 



3.2 Michel's monopole solution 

Before considering pulsar magnetospheres with a CKF-type 
structure, we check whether our routine works well. As a 
test, we try to get an exact solution of the pulsar equation 
call ed the Michel's monopo le solution which is given by (see, 
e.g, lOgura fc Koumall2003f ) 



*(R,z) = C I 



VR 2 + z- 



47rJ P (\E') = -n F * ( 2 
St(R,z) = 0, 



(30) 

(31) 
(32) 
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Figure 1. (a) Contour plot of ty. The nearest solid line to the axis is the magnetic field line with ty = 0.05C, and then the other solid lines 
are drawn from there by adding ty = 0.05C. The magnetic field lines pass through the light cylinder, R = -Rlo smoothly, (b) Contour 
plot of the relative error between the numerical result and the Michel's monopole solution. Sty is defined as Sty := K^n — ^mV^m 
where ty^ an d * m are the magnetic flux function of our numerical result and the Michel's monopole solution, respectively. We plot 
the contours that are Sty = 10~ , 10 -2 , 10 — 3 , and 10 — 4 . The contour line close to the origin represents Sty = 10 — 1 and the relative 
error is decreasing with distance. Since there is the line with Sty = 10 -3 near the ^-boundary, the relative error is not monotonically 
decreasing but the order is about 10 — 3 at most. Thus, we see that the relative error is quite small except near the origin. Since the 
Michel's monopole solution diverge at the origin, the relative error near the origin is larger than the other region but it is about 40% at 
most. 



where C is a constant. The Michel's monopole solution rep- 
resents a rotating magnetic monopole. We perform calcula- 
tions in the finite domain (0,0) < (R,z) ^ (27?lc, 2-Rlc) 
divided by 80 x 80 numerical meshes. In addition, we im- 
pose the boundary conditions as — at the axis R = 0, 
9/C = 1 at the equator 2 = 0, and ROr^ + zd z ty = 
at the outer boundaries. Then, we give the poloidal current 
function as Eq. (|31jl . Furthermore, we give rj and a in STnew 
as ?7 = 50 and a = 0.1 for this calculation and start our rou- 
tine with a non-zero initial toroidal current density. Results 
are shown in Fig. [T] We see that our result agrees with the 
Michel's monopole solution and it is smooth over the light 
cylinder. 



3.3 CKF-type solutions 

We study pulsar magnetosphere models that have a CKF- 
type structure. Unlike the CKF method, we need to deter- 
mine the functional form of ip (<]/). In order that, by in- 
specti ng the obtained polo id al cur re nts inlContopoulos et al.l 
1999h: " 



lOgura fc Koiimal (|2003h: iGoodwin et al.l (|2004) ; 
iGruzinovl (|2005h ; iTimokhinl (|2006h . we model the poloidal 
current Ip(ty) by using a cubic function. We describe our 



numerical model and show results. 



3.3.1 Boundary Conditions 

In this work, we consider a finite rectangular domain < 
R ^ fimax and < z ^ z ma x to solve E g. (1241) and impo se 
the boundary conditions as discussed in ITimokhinl (|2006l ). 

Firstly, let us consider the z-axis, R — 0. Since magnetic 
fields at the axis should not diverge, there is no magnetic 
flux at the axis. Then, the appropriate boundary condition 
is given by 



Then, at the equator, z = 0, in order to obtain a mag- 
netosphere with both closed and open magnetic field lines, 
we give the following boundary condition: 



d z V(R,0) 
9(R,0) 



for < R < i?LC, 



T„ 



for i? LC < R < Rn 



(34) 
(35) 



where f op is a constant. It is expected that the magnetic 
fields are closed for ty > ty op . Therefor, f op gives the value 
of the last open field line. 

Then, near the origin, we expect magnetospheres to be 
dipole which is given by 



mR 2 



(R 2 + 



for (0,0) < (R,z) < (Rs,z s ), 



(36) 

where m is a magnetic dipole moment and 7?s and zs rep- 
resent the surface of the star. We give this dipole's value at 
the surface and the interior of a neutron star. 

Lastly, at the outer boundaries, we impose boundary 
conditi ons with the mon opole magnetic field line as consid- 
ered by ITimokhinl I 



Rd R ty + zd z * = 0, for R = R n 
and < R ^ J? m „, z = z max . 



< z < z n 



(37) 



#(0,z) = for < z ^ Zn 



(33) 



For the CKF solution, other boundary conditions a re im- 
posed at the outer boundaries (|Ogura &c Koii ma 2003), how- 
ever, thei r results are quit e similar to the CKF solution ob- 
tained in ITimokhinl (|2006h in a region over the light cylin- 
der. It is expected that possible solutions are not sensitive 
to outer boundary conditions at least near the neutron star. 



3.3.2 Model of the Poloidal Current Jp(*) 

In our method, we fix the functional form of 7p(4'). In the 
previous works for the CKF solution, the obtained poloidal 
currents are similar. By inspecting their poloidal currents, 
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we assume that Iplp is given by 

16n 2 I P Ip = 16tt 2 A 2 *(* - 9 Ict ){9 - *o P ), 



(38) 



for ^ 9 < 9 o P where A and "frct^ ^o P ) are constants. 
Moreover, for 9 ^ f op , we assume that Iplp = because 
there is no poloidal current in the closed field zone. Since 
Iplp — (Jp)'/2, we can easily integrate Eq. ()38|) for ^ 



9 ^ ^o P and obtain 



47T/p(*) = ±^9^9* _ a^rct + ^op) ^ + 

(39) 

for ^ ^ < 9 op where we imposed the boundary condi- 
tion as lp(0) = so as to make the poloidal current den- 
sity regular at the axis. For 9 ^ ^op, which would be the 
closed field zone, we impose Ip(9) = 0. In this model, there 
are three parameters, A, 9 Ict , and 9 op . ^ret should satisfy 



^op/2 ^ 9 re t ^ ^ op not to make Eq. (|39|) be imaginary for 
«S * «S *op. 

'I' op represents the last open field line. Moreover, as- 
suming that the pulsed emission of pulsar originates from 
the polar cap, this parameter determines the emission area 
in the polar region near the neutron star surface. Let 9 P be 
the polar angle of the last open field line on the star sur- 
face. Since we put the dipole configuration as the boundary 
condition on the star surface, 9 P is given by 

. 1 „ *o 



sin f D = 



-°p 
9 S 



(40) 



where = m/Rs. When a rotating dipole is considered, 
^op is taken as ^ op = m/R^c =: that corresponds to 
the value of the magnetic flux on the light cylinder at the 
equator. Since Rs <C Rlc, &p is approximately given by 
#p ~ \J Rs /Rlc =■ do- Thus, in the case of our model, we 
find P ~ V*op/*o 9 . 

Let us discuss the meanings of the other parameters. 
From Eq. ((8|, the toroidal component of the magnetic field, 
Bt, is given by 



_E>t = 



2Jp(^) 
R ' 



(41) 



Therefore, roughly speaking, 7p expresses the toroidal mag- 
netic field. Thus, the parameter A determines the strength 
of the toroidal magnetic field; the choice of sign determines 
the direction of the toroidal magnetic field. Since Ip appears 
in the form of Iplp in the pulsar equation, the choice of its 
sign is not essential. We, hence, choose the minus as other 
authors did. 

Then, vE^t represents the existence of return currents. 
From Eq. (|38[) . we can easily see that 7p7p(^ r ct) = 0, and, 
from Eq. (|39[) . Ip(9 le t) 7^ in general. Namely, Ip becomes 
zero at 9 — 9 Tet , that is, the sign of Ip changes at 9 = 9 Tet . 
From Eq. (|12[) . the poloidal current density Jp is given by 



Jp = IpBp 



(42) 



Thus, since Jp is proportional to Ip, the return currents flow 
for * rot < * < * op . 

Thanks to fixing the functional form of the poloidal 
current, we can perform calculations with different types 
of the poloidal current with fixed boundary conditions. In 
addition, we can control the existence of a current sheet. 
Choosing ^ re t as 9 Tat = 9 op /2, it is easy to check that 
Ip(9 op ) = Oin this case. Therefore, we do not need a current 



sheet as introduced in the CKF solution. Moreover, tuning 
^rct = ^ op, we can consider the case that there exists only 
a current sheet, that is, there is no return current except the 
current sheet. 



3.3.3 How to give the parameters 

We model the poloidal current function Ip{9) by using a cu- 
bic function and it contains three parameters A, 9 Te t, and 
9 op . Although we try various sets of these parameters, the 
iteration does not always converge, which is expected by 
the regularity condition of the pulsar equation at the light 
cylinder. Therefore, we should give suitable parameters for 
calculations. Here, we discuss how we determine the param- 
eters A, 9 rc t, and ^ p- 

Let us introduce the ratio r as 



(43) 



Inter e stingly, in the previous works ( Contopoulos et all 
Il999l : lOeura fc KoiimJliooH : iTimokhinl [20061 ). the ratio is 
given by almost the same value that is about 0.8 though ^op 
is different from each other 0. It is expected that the ratio 
r is an essential parameter for the CKF solution. Therefore, 
we give the ratio r rather than 'I'ret- Since ^ re t should be in 
the range 9 op /2 ^ 9 Tet ^ ^o P , the range of the ratio r is 
0.5 < r < 1. 

Next, let us consider the parameter A. We investigate 
the poloidal current function of the CKF solution and find 
that it is quite similar to the Michel's monopole one's near 
9 = 0. The Michel's monopole solution gives that 7 P (0) 2 = 
f2 P /47r 2 , and therefore we adjust our model to the Michel's 
monopole solution near 9 = by using A as 



(44) 



Finally, giving ^ op and r, 9 Tet and A are automatically given 
by Eqs. (|43|l and (|44[1 , Examples of the poloidal current func- 
tion are shown in Fig. [2] We fix the ratio r in our numerical 
routine, and ^ op is determined not only to make the iter- 
ation converge but also to avoid unphysical solutions. For 
example, some solutions have closed field lines beyond the 
light cylinder, which would be unphysical. Our method can 
not always give a solution which is physically reasonable 
even if the iteration converge. We determine 9 op to avoid 
such unphysical solution. 



3.3.4 Results 

We want to construct a pulsar magnetosphere model with- 
out a current sheet which appears in the CKF solution. The 
existence of the current sheet and the amount of the return 
current can control by the ratio r in Ip(9). We perform cal- 
culations for various r in the finite domain (0, 0) < (7?, z) ^ 
(27?lc, 27?lc) divided by 80 x 80 numerical meshes. For the 
calculations, we give the value of r\ and a in St as r\ = 50 



1 The ratio 



Contopoulos et al.. 

Ogura fc Koiimal 1200$) , and 0.80 for 9 
d2006h . 



abou t 0.79 for \P p = 1.361-0 in 
|l999h , 0.82 for * op = 1 .669 in 
1.23*o in lTimokhinl 
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Figure 2. Poloidal current distribution of our model. The left panel and the right panel show 16n 2 IpI p and —Anlp for the ratio r = 0.5, 
0.6, 0.8, and 1.0, respectively. The case that r = 0.5 represents no current sheet model. Then, the case that r = 0.8 is the model like the 
CKF solution. The case that r = 1.0 is no return current model except a current sheet. 




R/Ri 



LC 



R/R 



LC 




R/R, 



LC 



R/R, 



LC 



Figure 3. Contour plots of the magnetic flux *. We show the poloidal magnetic field structures for r = 0.5, 0.6, 0.8, and 1.0. The thin 
dashed and solid lines in each panel represent the magnetic field lines. The nearest thin dashed line to the axis is the magnetic field line 
with if = 0.05'I'o! arltl then the values of the thin dashed lines and the thin solid lines are increasing from there with 0.05*1/0 an d 0.5*0i 



respectively. The thick solid line in each panel represents the last open field line: = 1.225*o f° r r = 0.5, 



1.234*o for r = 0.6, 



1.236*o for 



0.8, and * Q 



1. 234*0 f° r r = 1-0. The return current flows between the thick dotted line and the last open 



field line in each panel except the case that r = 1.0. In the case that r = 1.0, there is no return current except the current sheet. 



and a = 0.1. In addition, we assume that the star is repre- 
sented by the region of (0, 0) < (R,z) «S (0.057?LC,0.05i?Lc), 
which means that the pulsar period is about 4 msec assuming 
the star radius to be 10 km. We perform our numerical rou- 
tine with the convergence criterion discussed in Appendix iBl 
Results are shown in Fig. [3] from which we see that the 



magnetic field configuration near the neutron star is similar 
to the CKF solution for all ratios. Meanwhile, their global 
structures are different. In the case that r = 0.8, i.e., the 
model like the CKF solution, its structure is similar to the 
CKF solution in the whole domain. On the other hand, in 
the case without a current sheet, i.e., r — 0.5, the magnetic 
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Figure 4. The R-dependence of the magnitudes of the poloidal and toroidal current density. The top panels represent those at z = 
0.25i?LCi z — 0.5-Rlc> and z = I.ORlc m the case without a current sheet, that is, r = 0.5 . Then, the bottom panels show the case 
that r = 0.8 which is the case like the CKF solution and with a current sheet. Because we show the magnitude of the electric current, 
there are some sharp points in the figures, which means the change of sign of the electric current. Then, the rapid declines of the poloidal 
currents in the left pcncls indicate the existence of the dead zone. 



field becomes jet-like near the z-outer boun dary inside the 
light c ylinder. This is similar to the result in lLovelace et all 
{2006) in which free boundary conditions are imposed at the 
outer boundaries that is different from ours. Their solution 
shows a disk wind near the equatorial plane in addition to a 
jet flow along the axis, but our result has a quasi-spherical 
structure outside the light cylinder rather than a disk wind. 

One might think that a large toroidal current density 
appears near the light cylinder in the case without a current 
sheet, but it seems to have no clear relation. We compare the 
magnitudes of the poloidal and the toroidal current density 
for r = 0.5 and 0.8. Results are shown in Fig. [4] We see that 
there is no large difference between the poloidal and the 
toroidal current density in the both ratios. From the case 
without a current sheet in Fig. [3] it is easy to guess that 
the z-component of the magnetic field is superior to the R- 
component of that in the region where the jet-like structure 
appears. This structure is similar to an uniform magnetic 
field which is an exact solution of Ampere's law (|24|1 in the 
vacuum the same as the monopole solution (|30p and whose 
magnetic flux fucntion is given by ^ oc R 2 . Thus, a large 
toroidal current density would not be needed to make a jet- 
like strucure. From Fig. 2] the toroidal current density does 
not seem large, and therefore the appearance of a jet-like 
structure would not be odd as a solution of Ampere's law 
(|24p though we do not answer how the jet-like structure is 
formed. 



4 PHYSICAL PROPERTIES OF THE RESULTS 

In this section, we investigate physical properties of our re- 
sults. Firstly, we calculate the energy loss rate which is de- 
termined by integrating Ip(ty). In addition, we compare the 
total energy of magnetosphere between models we obtained 



in this study. Secondly, we study the difference of the global 
magnetic field structure between the case with and without 
a current sheet. Lastly, we check if the force- free approxima- 
tion is satisfied or broken down for our results. 

4.1 Energy loss rate of the magnetosphere 

We obtain pulsar magnetosphsere models both with and 
without a current sheet. Here, let us compare the energy 
loss rates that determines the power of a radio pulsar. In 
the force-free system, the energy carries away in the form of 
the Poynting flux. The Poynting flux, P, is given by 

P = —E x B. (45) 

4-7T 

Since the Poynting flux satisfies divergence-free condition, 
the total energy carried away is given by a surface integral 
on a closed 2-surface A: 

W= ( P-dS. (46) 

J A 

After some calculations (see e.g. iTimokhinl (|2006l )). it can 
be written in the following form 

/■*o P 1 

W = 2 — fi F /p(#)eZ*, (47) 

Jo 27r 

where the factor of 2 came from the fact that there are 
northern and southern hemispheres Q. Finally, the energy 
loss rate is determined by integrating /p(^I'). From Fig. [2] it 
is easy to find that the area bounded by Jp(^I') is the smallest 

2 Our for mula of energy lo ss rate is slightly different from the 
formula in ITimokhinl feooeft . It comes from the difference of the 
definitions of ip and \P. It is just notation and there is no essential 
difference. 
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Table 1. The energy loss rate for our results. Wo that appears 
in the table is equal to — m 2 /4tt 2 R^q. In the case of usual dipolc 
radiation, the maximum power represents (2/3)Wo- 



Ratio 



: 0.5 



0.6 



0.8 



1.0 



W/W 0.5002 0.6797 0.8228 0.8909 



1.014 



0.2R LC <(R,z)<1.0R LC - 
1.012 °-2H uc <M<2.0R LC ; 



1.01 
1.008 - 
1.006 
1.004 - 
1.002 - x 

1 : 



0.5 



0.6 



0.7 0.8 
ratio 



0.9 1 



Figure 5. Comparing the total energies. We normarize each total 
energy by the case that r = 0.5. We estimate the total energies in 
the domains (0.2_R LC , 0-2R LC ) < (R, z) < (I.ORlc , LORlc) and 
(0.2ii LCl 0.2.R LC ) < (R,z) < (2.0fl LC ,2.0fl LC ). The difference 
of the total energy is about 1.2% at most and becomes smaller as 
wc take the integration region widely. 



for r = 0.5 that is the case without a current sheet, and then 
it monotonically increase with increasing ratio r up to unity 
if "i/op for all ratios are the same. For our numerical results, 
the values of ^ op are almost the same. Thus, the energy 
loss rate in the case without a current sheet would be the 
smallest in our model. In fact, we calculate the energy loss 
rates for our results and show them in Table [T] Considering 
pulsar physics, the energy loss rate affects the characteristic 
age of a pulsar. Our result implies the characteristic age of 
a pulsar depends on the existence of a current sheet and the 
difference is about the factor of 2 at most. 

In addition to the energy loss rate, we can also calculate 
the total energy of electromagnetic field which is given by 
the volume integration as 



B ) dV. 



(48) 



Results are shown in Fig. [5] Although the case without a 
current sheet has the least energy, the difference from other 
models is quite small. This is because the magnetic field 
structures at the inner and outer boundary are irrespective 
of the ratio parameter r in ip(\l/). The integrand is domi- 
nated by these boundary values. At the inner boundary, the 
energy density is much larger than outer part because we 
impose the boundary condition as the magnetic field struc- 
ture to be dipole, while at the outer boundary the volume 
element is much larger than the inner part: \B\ ~ 1/r 3 near 
the origin where r — \/B? + z 2 ; \B\ ~ 1/r at a distance 
due to the toroidal magnetic filed; and the volume element 
is proportional to r 2 . Therefore, the integrated energies be- 



come similar for all models. Although the difference of the 
total energies is quite small, the energy loss rate without a 
current sheet is clearly smaller than those with a current 
sheet. Thus, the model without a current sheet seems most 
plausible in our models. 



4.2 Three dimensional magnetic field structure 

The magnetic field structure is important when one consid- 
ers a mechanism of radiation from pulsar, e.g., curvature ra- 
diation. Therefore, it is worth showing the three dimensional 
structure of the magnetospheres. One magnetic field line 
can be drawn by integrating B which is given by Eq. (|15[) . 
Here, we show the three dimensional magnetic field struc- 
ture with two methods; the structure on z — const, plane, 
and the structure in a three dimensional box. In order that, 
we change the cylindrical coordinates (R, ip, z) to the Carte- 
sian coordinates (x,y,z) and integrate B on a (x, j/)-plane 
or in a (x, y, z)-box. We draw the magnetic field structure in 
the case that r = 0.5 and 0.8 and results are shown in Figs[S] 
and [7] We see that the magnetic field lines in the case with a 
current sheet are more twisted than those of the case with- 
out a current sheet. This is consistent with the discussion of 
the energy loss rate because the radiation gets stronger for 
the magnetic fields with smaller curvature radius. Because 
of the difference of Ip (*]/), the three dimensional structure 
is quite different from each other. 



4.3 Violation of the force-free approximation 

In our routine, the convergency is checked whether the 
change of the toroidal current function St is in a small 
residual. If St would hardly change, it would imply that 
the force-free condition is sufficiently satisfied because we 
give the next St by using the force- free condition (|25p . In 
the other word, the pulsar equation is satisfied. This is true 
inside and outside the light cylinder, but there is no guaran- 
tee at the light cylinder because we used the light cylinder 
condition (|21[) to give Stlc though it is not satisfied during 
the iteration. We expect that the light cylinder condition 
would be automatically satisfied when our routine is fin- 
ished. In fact, we can obtain the Michel's monopole solution 
as shown in §3.2\ However, whether the light cylinder con- 
dition is satisfied after the iteration converge is not always 
guaranteed in our method unfortunately. 

We check whether the pulsar equation is satisfied for 
our results and find that it breaks to some extent around 
the light cylinder. The violation of the force-free condition 
at the light cylinder in the case that r = 0.5 is shown in 
Fig. |SJ for example. Finally, we find that our results do not 
satisfy the pulsar equation near the light cylinder that is 
about 0.9-Rlc < R < I.IRlo- This violation is seen in 
all ratios in almost the same region, but the violation is 
most suppressed in the case r = 1.0. In addition, we tried 
other outer boundary conditions given bv lOgura fc Koiimal 
(|2003h and perform the same calculation. Although the vi- 
olation of the force-free condition is also appeared in this 
case, the degree of the violation is improved comparing with 
the boundary condition discussed in H3.3.1I Moreover, the 
width of the violating region seems to come from the pa- 
rameter a in Si-new which determines the variance of the 
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Figure 6. Magnetic field structure on a (x, y)-planc. The top panels represent the magnetic field structure in the case that r = 0.8, and 
then the bottom panels represent those of the case that r = 0.5. These figures show the magnetic field structure on the (x, y)-plane with 
Z = 0.25RlC; 2 = 0-5-Rlc> and z = 1.0-Rlo The circle m the center of each panel indicates the light cylinder. We see that the magnetic 
field structure with/without a current sheet is quit different from each other near the equator. The magnetic field lines in the case with 
a current sheet are more curved than those of the case without a current sheet. 




Figure 7. Magnetic field structure in a (x, y, z)-box. The top panels represent the magnetic field structure in the case that r = 0.8, 
and then the bottom panels represent those of the case that r = 0.5. Three straight lines and the black ball in the center in each panel 
indicate the x, y, and z axis and a neutron star, respectively. The magnetic field lines in the left panels arc drawn from various points: 
(R,<p,z) = (0.02iJ LC ,'t7r/i2,0.14i?Lc), (0.045iJ L c , "r/ 24 . 0.115R L c), and (0.06iJ L C , "r/ 24 . 0.12_R LC ) where i = 0,1,2,--- ,24. For the 
right panels, we give the initial points at (R,ip,z) = (0.25 sin dR-^Q, 0, 0.25 cos 0Rlc) where 9 is the angle from z-axis and it takes as 
= Mr/48 for i = 0,1,2, • •• ,24. 
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Figure 8. Violation of the force-free condition in the case without 
a current sheet. We estimate the magnitude of the pulsar equation 
II2UI at the light cylinder, that is, the light cylinder condition (3TJ. 
The vertical axis shows the magnitude along the light cylinder 
which is normalized by the summation of the absolute value of 
each term in ( 121)1 ). The dots are the data of our results and they 
are connected with line. 



Gaussian function. Giving other a, the width of the viola- 
tion is changed. To sum up, we could make the violation of 
the force-free condition improved by adjusting outer bound- 
ary conditions and the poloidal current function though we 
do not know how to give them suitably. 

Concerning to the astrophysical application, it is widely 
known that the fast-magnetosonic point is well-inside of the 
light cylinder for most pulsars other than the Crab pul- 
sar whose fast-magnetosonic point is nearly equal to or 
slightly further than the light cylinder. For example, the 
Lorentz factor of the fast-magnetosonic point can be writ- 
ten in jf m = ^fa where a represents the ratio of the Poynt- 
ing energy flux to the matter energy flux and is given by 
727 where n is the particle number density 
in the pulsar rest frame, m e is the ele ctron mass, and 7 i s 
the Lorentz factor of the bulk velocity (|Kirk fe Duffy|[l999l) . 
By eval uating the number density by the Goldreic h- Julian 
density (|Goldreich fc Julian|[l969l ; iKirk et alJl200Sh . we can 
estimate the fast-magnetosonic velocity at the light cylinder 
for various pulsars. In the case of the Crab pulsar, the ratio 
of Lorentz factor of fast-magnetosonic velocity and obser- 
vationally estimated value at the light cylinder is ~ 0.98; 
for other pulsars whose parameters are observed or well es- 
timated (|Tanaka fc Takaharal l20LlT ). the ratio is typically 
order of 1CP 6 , which means the fast-magnetosonic point is 
at the well inside of the light cylinder for most observed pul- 
sars. The force-free approximation would break at the fast- 
magnetosonic point so that the breakdown of our results at 
some point near the light cylinder is not serious problems 
and can be applied to many astrophysical pulsar magneto- 
spheres. In particular, in the case without a current sheet, 
our result has a jet-like structure inside the light cylinder 
and it gives us another possibility of CKF-type solutions. 



5 SUMMARY AND DISCUSSION 

We numerically studied stationary pulsar magnetosphereres 
in the force-free system both with and without a poloidal 
current sheet which appears in the CKF solution. In this 
study, we expressed the existence of the current sheet by the 
ratio r that represents the amount of the return current and 
constructed pulsar mangetosphere models for various ratios. 
In the case with a current sheet, in particular near r = 1.0, 
the structures on (R, z)-plane are similar to the CKF so- 
lution. Meanwhile, in the case without a current sheet, i.e. 
r = 0.5, it has a jet-like structure inside the light cylinder 
and becomes quasi-spherical at a distance, which is differ- 
ent from the CKF s olution but similar to the solution in 
I Lovelace" et all (|2006r i. Then, we studied physical properties 
of our results and showed that the existence of the current 
sheet affects the properties of pulsar magnetoshpere. First, 
we calculated the the energy loss rates for our results and 
showed that it depends on the ratio r. It would affect the 
characteristic age of a pulsar. Moreover, the model without 
a current sheet has the minimum energy and energy loss 
rate though there is no large difference for the total ener- 
gies. It sounds that the model without a current sheet is 
plausible. Next, we showed their three dimensional struc- 
tures both with and without a current sheet. From Figs. [6] 
and we can find that how curved the magnetic field lines 
hardly depend on the value of z in the case with a current 
sheet. In the case without a current sheet, the magnetic field 
lines are more bent as z increases. 

The difference of the magnetic field structure between 
the case with and without a current sheet suggests that the 
curvature radiation from the pulsar depends on whether a 
current sheet exists. In the case with a current sheet, it seems 
that the curvature radiation from the magnetosphere would 
uniformly occur in the whole region. On the other hand, in 
the case without a current sheet, occurring the curvature 
radiation would concentrate in a polar region. We show the 
angle between the magnetic field line and R or z axis in 
Fig. [5] These angles, <f) an< A Xi are given by 



cos — 



cos x 



B-e R 
\B\ ' 
B-e z 

\B\ ' 



(49) 
(50) 



where Sr and e z are a unit vector tangent to the R and 
z direction, respectively. The data are useful to calculate 
the curvature radiation, which i s a plausible mechanism 
for ob served pulsar radio emission (|Ruderman fc Sutherland! 

Eazsl). 

Our models do not satisfy the force-free condition near 
the light cylinder, however, concerning realistic astrophys- 
ical situation, the force-free approximation would break at 
the fast-magnetosonic point which is inside the light cylin- 
der. Therefore, our models could be used enough in realis- 
tic astrophysical situation. Moreover, the physical proper- 
ties discussed in this paper are the same even if we discuss 
the properties only inside the light cylinder. It is needless 
to say that effect of inertia of plasma is important in the 
wind zone and dissipation should also be concerned. There- 
fore, another treatment beyond the force-free approxima- 
tion should be concerned, which is our future work. It is 
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Figure 9. The angle between the magnetic field line and R or z. The top and bottom panels show the case that r = 0.8 and r = 0.5, 
respectively. We calculate the angle cos</> and cosx with z = 0.25-RlCi 0.5i?LC> an< l 1-0-Rlc f° r tne both ratios. 



an open and interesting question whether the jet-like struc- 
ture which is seen in the case without a current sheet ac- 
tually appears when one concerns the effect of inertia. In 
the framework with stationarity, we can not answer if a 
magnetoshpere without a current sheet appears and if the 
j et-like structure is re alized. The Chandra X-ray satellite 
ijWeisskopf et al.ll2000T ) has recently found the jet-like struc- 
ture at the polar region of the Crab pulsar wind nebula. Re- 
cent relativistic magnetohydrodynamical studies have shown 
that this jet can be resulted from compres sion by the mag- 
netic hoop stresses or large scale vortexes jDel Zanna et al.l 
12004 iKomissarov fe Lvubarskvll2003l ; lLvubarsky||2002l ). Al- 
though it is difficult to show that our jet-like structure is 
responsible for the Crab jet due to the force-free approxi- 
mation, we think that our results can be a new probable 
initial condition for the RMHD studies of pulsar wind neb- 
ula, instead of the split monopole model. 
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APPENDIX A: DERIVATION OF THE NEW 
TRIAL TOROIDAL CURRENT 

In this appendix, we will derive Sinew First of all, let us 
derive Sr-in, STout, and Stlc- From Eqs. p4l) and ([25]), we 
have 

R 2 nlS T (R, z) - 2mld R <H + 16tt 2 /p/p = S T (R, z), (Al) 

St(R,z) in the right or left hand side of the above equa- 
tion can be regarded as a trail toroidal current Serial z). 
Thus, we can make the following two equations 0: 

S Ti n{R, z) = i {(1 + i? 2 fi|)S Ttr ial(ii, Z) 

-2fi|i?c> H * + 167r 2 /pJp}, (A2) 
S To ut{R,z) = (1 + R 2 ^)- 1 {25 Ttr iai(i?,z) 

+2nlRd R V - 167r 2 / P Jp} , (A3) 

5*Tin or STout can be regarded as the new toroidal current in 
our method. Unfortunately, the iteration only converges if 
the numerical domain is inside or outside the light cylinder. 
Then, in order to use our method in a domain including 
the light cylinder, by using Stiu, STout and the hyperbolic 
tangent, we make the following function STsidc z): 



STside(R, Z) 

_ 1 + tanh^D) 



S Tin (R,z) + 



tanh(?7_D) 



Sto 



t(R,z). 
(A4) 



Thanks to the hyperbolic tangent, SWae almost becomes 
Stiii and STout inside and the outside the light cylinder, 
which is expected to make our iteration converge across the 
light cylinder. 

STsidc has useful form to make the iteration converge, 
however, we still need modification. From Eqs. (|A2|) and 
(|A3[I . it is easy to see that both of STin and STout be- 
come STtriai at the light cylinder, that is, STsidc (-Rlc, = 
STtriai(-RLC, z) . Thus, it forces the value of STside at the 
light cylinder to fix during the iteration. We do not know 
the exact value of the toroidal current at the light cylinder 
a priori, and therefore we need some modification to STsidc 
in order to change the value at the light cylinder during the 
iteration. 



3 To avoid the divergence at the axis R = 0, we added St both 
sides of Eq. jAll l before making Eqs. \A2l and J A3ti ■ 



Let us return to Eq. (|A1[> . The exact form of the toroidal 
current is given by 



S T (R,z) 



-gnigggg + i67r 2 jp/p 
(i-i? 2 n|) 



N 

15 



(A5) 



Let the regularity condition at the light cylinder, N = at 
R = Rlc, be satisfied. Then, by using l'Hopital's rule, the 
toroidal current at the light cylinder is given by 



N 
T) 



d R N 



9 R D R=R LC 
R=Rlc 



-RZc^v 2 (^ 2 IpIp)'9r^\r=r l 



(A6) 

This is Stlc(-Rlc, z) and gives the exact value at the light 
cylinder formally. 

Finally, by using the Gauss function, we have the fol- 
lowing form of the new trial toroidal current: 

(R,z) 

1 + tanh(,D) ST;n(7?) z) + l-tanh(^) gTout( ^ g) 



2 

x ll - c 



1 1 — D 2 

| + Stlc(-Rlc, z)e 



(A7) 



Thanks to the Gauss function, JTnow becomes Jtlc at the 
light cylinder. 



APPENDIX B: CRITERION OF THE 
CONVERGENCY 

We discuss the criterion of the convergency of our calcula- 
tion here. During the iteration, we solve Ampere's law (|24|l 
with various St that is an elliptic type differential equation. 
Solving Eq. (|24jl , we use the conjugate gradient method and 
the convergecy obeys the method. In addition, we should 
check whether St converges and it is determined as follows. 
Let us introduce the following quantity 



Hi 



+ d 2 R * N -±-d R * N + d 2 * N , (Bi) 

R 



where N indicates the number of iteration times, and S^ +1 
is determined from Eq. (|26[) with . Then, we integrate 
the magnitude of Hp over the numerical domain, A, as 



\Hp\dS = H, 



(B2) 



where S is the area of A. We use H to judge whether St con- 
verges. If H can be numerically ignored, there is no difference 
between S* +1 and St • In our calculations, the derivatives 
of ^ are written with the second-order accuracy. Therefore, 
we continue our routine until H can be ignored with the 
second-order accuracy. We show an example in Fig. IBll and 
see that our method that gives a next toroidal current den- 
sity as Eq. (|26p shows the second order convergency at least. 
In this calculation, we give rj and a in Stiicw as a = 4Ah 
and T) = 0.5ct -2 where Ah(— AR — Az) is the grid interval. 
It could be expected that Eqs. (|24p and (|25p are numerically 
satisfied simultaneously with ^ N and after St converges. 
However, it is not guarantee at the light cylinder because we 
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Figure Bl. The plots of H in the case that the outer boundaries 
are set at -R max = 2max = 2_Rlc- The vertical and the horizontal 
axis are H and Ah 2 , respectively. The criterion of the convergency 
we imposed are satisfied with the second-order accuracy. 

do not impose the light cylinder condition during the iter- 
ation. Although whether the force-free condition is satisfied 
at light cylinder is not guarantee in our method, our method 
gives a convergence solution of Ampere's law with St giving 
by (H)}. 

This paper has been typeset from a TfrjX/ J^Tj^X file prepared 
by the author. 
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